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. Abstract Hybrid models of chemotaxis combine agent-based models of cells 

with partial differential equation models of extracellular chemical signals. In 
this paper, travelling wave properties of hybrid models of bacterial chemo- 
taxis are investigated. Bacteria are modelled using an agent-based (individual- 
based) approach with internal dynamics describing signal transduction. In ad- 
' dition to the chemotactic behaviour of the bacteria, the individual-based model 

also includes cell proliferation and death. Cells consume the extracellular nu- 
trient field (chemoattractant) which is modelled using a partial differential 
equation. Mesoscopic and macroscopic equations representing the behaviour 
of the hybrid model are derived and the existence of travelling wave solutions 
. for these models is established. It is shown that cell proliferation is neces- 

^ ' sary for the existence of non-transient (stationary) travelling waves in hybrid 

models. Additionally, a numerical comparison between the wave speeds of the 
i-^ ' continuum models and the hybrid models shows good agreement in the case of 

weak chemotaxis and qualitative agreement for the strong chemotaxis case. In 
the case of slow cell adaptation, we detect oscillating behaviour of the wave, 
which cannot be explained by mean-field approximations. 

> 
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1 Introduction 

The wavelike spread of cell populations plays a fundamental role in many bio- 
logical processes, including development [22, wound healing [38] and tumour 
invasion |16] . Bacterial populations show similar phenomena, with the pio- 
neering studies of Adlcr [I] confirming the capacity of an E. coli population to 
form travelling bands via chemotaxis to extracellular signals. Mathematically, 
the extent to which chemotaxis can generate and sustain stationary travelling 
bands has motivated a number of studies, including the Kcllcr-Scgel model of 
Adler's experiments which is written in the form of coupled partial differential 
equations (PDEs) [50]. This early model necessitated a biologically unrealis- 
tic singularity in the chemotactic sensitivity to generate stationary travelling 
waves: a requirement that allows bacteria behind the wave to acquire infinite 
speeds and to avoid "dropping-out" , an effect that leads to gradual dispersal 
of the band [40lll5]. 

This singularity requirement can be circumvented by incorporating other 
processes. The well known Fisher's equation [M] demonstrates travelling waves 
in systems coupling diffusion with logistic growth terms |14] . Parabolic chemo- 
taxis models with non-singular sensitivities but incorporating either logistic 
[^[^[50] or non-logistic [^155] growth terms also admit travelling wave so- 
lutions. Other studies have shown that introduction of more complex nutrient 
terms can give rise to travelling waves, even when growth is absent [341135] . An 
experimental system which also included two chemicals - a chemoattractant 
and a nutrient source - was presented in [SJIT], with stationary or transient 
travelling waves obtained according to t he formulation of the model [51140) . 
Travelling waves in chemotaxis models have also been recently studied in [261 
[25] ; we also note the articles [19] and ]37] for a review and analysis of travelling 
waves in PDE-based models. A comparison between mesoscopic (hyperbolic) 
and macroscopic (parabolic) PDEs has been presented in [??] ■ 

Relatively little exploration has been conducted into travelling wave for- 
mation for chemotactic models extending beyond PDE systems, in particular 
those introducing terms to account for the inherent noise of biological systems. 
One exception is the study of [S] , in which a multiplicative noise term was in- 
troduced into the Keller-Scgcl model and the existence of travelling waves has 
been demonstrated within this setting. Hybrid models, in which an individual- 
based model for bacterial behaviour is coupled to a continuum description of 
extracellular signals, naturally introduce stochastic effects and will be the focus 
of the present paper. Such a hybrid model was formulated in |15j where it was 
shown that under finite cell speeds only transient travelling waves formed, even 
with singular chemotactic sensitivity. The individual-based model was formu- 
lated in terms of the velocity-jump model with internal dynamics [T2lll3ll4T] 
and, in this paper, we extend the model in [15] to incorporate proliferation 
and death of bacteria. We analyse this system numerically and analytically 
with respect to its travelling wave properties, employing the biologically in- 
spired chemotactic sensitivity presented in [40j and a linear growth term. We 
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show that stationary travelhng waves can be observed even in the absence of 
chemotaxis, although wave speeds are substantially increased in its presence. 

The organisation of the paper is as follows: the full hybrid model is pre- 
sented in Section [2] along with illustrative simulation results, while the corre- 
sponding continuum equations are derived under certain assumptions in Sec- 
tion 131 in Section U these continuum equations are analysed with respect to 
travelling wave properties; in Section [5] where a computational analysis and 
comparison of the models is presented; finally, we discuss our observations in 
Section [6l 



2 Hybrid model of bacterial chemotaxis 

In this section we formulate the hybrid model of bacterial chemotaxis which 
will be investigated in this paper. The model is motivated by the behaviour 
of the bacterium E.coli and, in its most general form, includes cell movement, 
sensing and response to a chemical signal, consumption of the chemoattrac- 
tant, cell proliferation and death. However, for analytical tractability, we will 
also explore simplified hybrid models which exclude some of these processes. 
Bacteria are modelled as agents with internal dynamics that represent the sig- 
nal processing and response of each individual while the extracellular chemical 
is modelled using a PDE to describe its spatio-temporal concentration. The 
mathematical framework and simulation techniques are reviewed in |15] . We 
consider the model in an effectively one-dimensional domain representing a 
long but narrow tube, similar to the experimental set up considered in [T]. 

The motion of E. coli bacteria is controlled through the coordinated rota- 
tion of flagella distributed over the cell surface [5] . Counterclockwise rotation 
generates a propulsive bundle that results in straight line motion of the bac- 
terium - a so-called "run" [5]. Alternatively, clockwise rotation results in the 
outward flaying of flagella and a "tumble" - rotation with insignificant dis- 
placement. At the end of each tumble the bacterium chooses a new direction 
of movement, seemingly at random, and returns to the run phase. The lengths 
of the individual phases are independent from each other and distributed ex- 
ponentially, yet they can be influenced by internal dynamics [2]. 

Internal dynamics of the E. coli bacteria possess two principal features [4] : 
a quick excitation phase followed by slower adaptation. Specifically, changes in 
the extracellular signal concentration lead to quick excitation of the internal 
metabolism, signified through altered chemical concentrations inside the cell. 
Following excitation the internal concentrations revert slowly to normal in an 
adaptation process, even when the external signal remains at the raised level. 

2.1 Velocity jump model with internal dynamics 

Run- and- tumble dynamics are aptly modelled as a velocity-jump process |31l 
[T^ . We denote by Na{t) the number of bacteria (agents) in the system at time 
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t. The current state of the i-th individual, i = 1, 2, . . . , Na{t), will be described 
using its position Xi € R, its velocity Vi = ±s G R and a set of internal state 
variables € M™ that represent the states of components in the intracellular 
signal transduction network. 

Here we concentrate on a cartoon version of the internal dynamics of bac- 
teria written in terms of two internal variables [32012] . i.e m = 2. Internal 
variables and y'^^'' are governed by the equations 



(2.1) 



dt te 

dy(2) S'(a-(t),t) -2/(2) 



dt ta 

where te is the excitation time, ta is the adaptation time, te <C ta and S{x{t), t) 
is the concentration of chemoattractant at the position of the bacterium x{t) 
at time t. Furthermore, bacteria move with the velocity Vi = ±s governed 
through a velocity jump process with a turning frequency A = A(y) that 
depends on the internal dynamics. In this paper, we will use the biologically 
motivated nonlinear turning kernel developed in [40] . Hence, the full model of 
one individual over (a small) time step At can be written as: 

x{t + At)=x{t)+v{t)At, (2.2) 
^ ^ r -t-(t), with probability \{y{t)) At, 
y i)(f), otherwise , 

X{y{t)) =xJl- , (2.4) 



«+|2/(i)(i)|; ' 

.^ (^^,^ S(x(t),t) - y(^'> (t) - y^^^t) ^ 
y'-^'it + At) ^ 2/(i)(t) + ^ ^ ^' ^ a — U. a — \-LAt, (2.5) 

2/(2) (i + ZiO = i'Kt) + ^t, (2.6) 



where Aq and k arc positive constants. 

In addition to the behaviour of an individual bacterium we define a signal- 
dependent proliferation function h{S) : R"*" i— >■ R. We thereby interpret a 
positive value of h(S) as a proliferation rate, meaning that in the infinitesimal 
interval [t, t + At) a bacterium at position x generates an exact copy of itself 
with probability h{S{x{t),t)) At. Similarly, a negative value of h{S) means 
that the bacterium disappears (dies) with the probability —h{S{x{t),t)) At. 
In this paper, we will use the following form for the proliferation rate h(S): 

h{S) = a{S - Se) , (2.7) 
where a and Sc are positive constants. 
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2.2 Evolution of the extracellular chemoattractant 

For the extracellular signal S{x, t) we formulate a PDE that incorporates dif- 
fusion (with diffusion constant Ds > 0) and signal consumption by bacteria, 
the latter with signal dependent rate k{S) : R+ R+. The equation for S 
therefore takes the form 

2 — 1 

For the remainder of the paper we employ a linear form for the consumption 
function k{S): 

k{S) = l3S, (2.9) 

where /? is a positive constant. 



2.3 Illustrative example 

The hybrid model framework presented in Sections 12.11 and 12.21 includes es- 
sential features of the more complicated hybrid chemotaxis models formulated 
in [T0ll39] . In this section we numerically show that these processes can give 
rise to travelling waves. For the numerical simulation we employ techniques 
described in |15] . In particular, for the extracellular signal S{x,t), this means 
that the simulation is performed on the one-dimensional domain [0, L] with 
initial condition S{x,0) = Sao > and zero-flux boundary conditions. We 
consider M + 1 regularly spaced grid points rj = j Ax, j = 0, . . . , M, where 
Ax = L/M and the values of S{xi,t) are advanced by a small time step At 
and a forward Euler update rule: 

S{r„t + At) = Sir„t) + Ds ^t ^ir^^^^t) + S{r,^,,t) - 2Sir„t) 

(Ax) 

N^it) (2.10) 
-k{S{r„t))At J2 K(r,~xm- 

In the above K : M. ^ IR+ is the symmetric, normalised and non-negative 
kernel 

e 



V 27rcr^ 



20-2 



where the kernel width cr is a positive real number. Here, K(rj —Xi) represents 
the influence a bacterium at position Xi has on grid point j. 

The simulation of the individual bacterium is given in the full system (|2.2p - 
(|2.6p and complemented by the birth and death processes described in Sec- 
tion where we use the same time step At as in (|2.10|1 . To calculate the 
necessary off-grid values of extracellular signal, we linearly interpolate from the 
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two nearest grid points. We further simplify the system (j2.2p - (|2.6p by exploit- 
ing the separate time scales for excitation and adaptation (i.e. te <C ta): specif- 
ically, wc assume the update equation (|2.5p for y^^^ is in a quasi-equilibrium, 
which is identical to the assumption te ~ 0. The value for 'ip-^ can therefore 
be calculated by 

V^^\t) = S{x(t)^t)-y^^\t). (2.11) 

Illustrative results are presented in Figure [T] For this simulation, Afa(O) = 10^ 
bacteria were initialised at positions a;i(0), randomly generated as the absolute 
value of a Gaussian random variable with variance much smaller than the 
domain length L. The initial velocity (direction of movement) is generated 
uniformly at random and initial values of the extracellular signal and internal 
variables are taken as 

yl'^(0)=5oo, 2/f^(0) = 0, fori = l,2,...,iVa(0), 
S'(a;, 0) = 5*00 forxe[0,L], 

where 6*00 = 1. We simulate the system until time Tgnai = 100 and plot both 
the distribution of bacteria and concentration of chemoattractant S in Figure 
[TJa). We also estimate the wave speed as a function of time in Figure [TJb). 
We clearly see formation of a travelling band of bacteria, moving rightwards 
with average speed v = 0.51 (plotted as the dashed line in Figure [ijb)). 

Influence of the growth term 

To investigate the influence of the growth term on the existence of travel- 
ling waves, we simulate the full hybrid model (|2.2|) - (|2.6|) and (j2.10|) including 
(a = 1) and excluding (a = 0) growth and death processes. Wc use identical 
parameters to those described above and present the results in Figure [5] In 
Figure [2^a) the position of the wave front (defined as the right-most position 
for which S{x) < 0.9) is compared. The full hybrid system (dashed line) gen- 
erates a straight line, indicating a wave moving with constant speed. While 
the system excluding growth and death (solid line) moves with a similar initial 
speed, speed is gradually lost over time: the shape of n{x, t) at different times 
for this case is shown in Figure [2Kb). We clearly see that no true travelling 
wave forms, with many agents being left far behind the wave front, leading to 
its slowing down. Thus, we can interpret growth and death terms in terms of 
a stabilising role on the wave profile: although not all agents can keep up with 
the wave, new agents are constantly created at the front and the agents that 
drop out eventually die, resulting in a travelling band of agents. 

3 Prom hybrid models to macroscopic PDEs 

In this section we derive macroscopic PDEs for the spatio-temporal density of 
bacteria n{x, t) at given position x € M and time t > 0. An implicit assumption 
of the derivation is spatial independence of bacteria, which allows formulation 
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Fig. 1 Numerical solutions of the hybrid chemotaxis model I l2.2l l-l l2.6l l and 112.101 1 and PDE 
System A ll3ll-ll331l. 

(a) Wave form for the hybrid model after time t = 100. Solid line: estimated density of 
bacteria, dashed line: extracellular chemical signal S, 

(b) Measured speed of travelling wave (solid line). Dashed line denotes the average speed. 

(c) Wave form for PDE system A after time t = 100. Solid line: estimated density of 
bacteria, dashed line: extracellular chemical signal S. 

(d) Measured speed of travelling wave (solid line) for PDE System A. Note that the spike 
near t = is a product of the wave speed calculation method. 

The dimensionless parameters are: a = /3 = s = l, Sc = 0.5, Soo = 1, = 10""^, 
Ax = 0.25, L = 100, Ao = 10, k = 0.01, Ds = 0, ta = 0.1, a = 0.5, 



of a continuous mesoscopic system. We then use results from [12] to obtain 
the macroscopic equations. To illustrate the successive formulation of models 
we construct two systems of PDEs - denoted System (A) and System (B) - 
to be referred to in the remainder of the paper. 



3.1 System (A) 

We define the mesoscopic densities p^(x, y^^^ , t) for left and right- moving bac- 
teria, depending on their position a; G M, their internal variable y'^' S M and 
t > 0. If the signal profile S = S{x,t) was uninfluenced by bacteria, densities 
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Fig. 2 Numerical solutions of the hybrid chemotaxis model 112. 2| I -1 12.61 1 and I I2.10I I without 
growth and death terms. 

(a) Comparison of position of wave front over time. Solid line: without growth/death (a = 
0), dashed line: with growth/death (a = I). 

(b) Wave form at different times during simulation with a = 0. From left to right: t = 
20,40,60,80. 

Remaining parameters as in Figure U\ 



would satisfy the following system of hyperbolic PDEs: 
dp+ , , d f S{x,t) - yf-^^ ^ 



dp- , d fSi.,t)^yi^[^ = Xp^-Xp-+KS{.,t))p-, 



dt dx 9?/(2) V ta 

(3.1) 

where A is defined in (|2.4p which, under (|2.1ip . can be simplified to 

The signal dynamics is described by (j2.8p which can be rewritten in terms of 
p^ as 

w = " jy ^p'^'^y^'^ ■ (3.3) 

We denote the system of equations (j3.ip - p.3p as System (A). 

The system p.ip (for the one-particle distribution) can be derived by inte- 
grating the probability distribution function p(xi , wi , yi ; X2 , W2 , y2 ; ■ • ■ \S{x,t)) 
for the many particle system, utilizing the fact that the movement of individ- 
uals are biased by the signal function S{x,t), but independent to each other. 
However, for the hybrid chemotaxis models described in Sections 12.11 and 12.21 
individual bacteria interact via the extracellular signal S which complicates 
the derivation of p.ip . In [11], a kinetic description has been derived for a 
model of interacting locusts, using a modified version of the BBGKY hier- 
archy from the classical kinetic theory of gases [H]. The system we consider 
here is much more complicated to analyse than the locust model studied in 
[11) . due to the variable number of bacteria and internal variables. Thus the 
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kinetic description p.ip can only be considered as an approximation to the 
one particle distributions of the interacting system. 

The capacity of the above mesoscopic system to generate travelling bands 
analogous to those observed in the hybrid model is illustrated in Figure [Hc)- 
(d). For details of the numerical method employed for this and other simula- 
tions of the continuous model, we refer to |40| . The qualitatively and quan- 
titatively close correspondence in solutions under equivalent parameters and 
initial conditions corroborates the use of the above approximation. 



3.2 System (B) 

We consider a macroscopic model in this section. Define the macroscopic den- 
sities 

p^ix,t)= f p±(x,y(2),t)dy(2), (3.4) 
Jr 

and let them satisfy the following system 



dp~ dp- , , [dS\ , fdS\ _ 



(3.5) 



dt dx \dx J V dx 

where the turning rates are given by 

A^ = AoflTx§^) with X- , nt°9X. V (3.6) 
V ox J KXa[l + 2Xota) 

Using p.4p . equation p.3p can be written as 

We will denote p.Sp and (|3.7p along with the definition of A* in (|3.6p as 
System (B). According to the analysis in [12[HT| . System (B) is quantita- 
tively consistent with System (A) when the external signal S{x) changes slow 
enough such that cells are close to their fully adapted state, in which case cell 
movement is only moderately modified by the signal. 



In the rest of the paper, we assume diffusion of extracellular signal to 
occur on a much slower time scale than the active motion of the bacteria, 
hence Ds = 0. The number of parameters of the above models can be reduced 
by setting s,Soo,ct,f3 to one through rescaling. We show this in detail for 
System (B) as follows. Rescaling the variables S = 55*00, = P^a5oo//3, 
t ~ i/(a5oo), X = xs/{aSoo) and the parameters Sc = ScS^o, Aq = XoaSoo, 
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taking (j2.7p and substituting into System (B) wc obtain, after dropping hats 
for notational simplicity, 

= ^*('Dp*-^-('Dp- + (S-S.)p-. (3-8) 





9p+ 


~dr 


ox 






dt 


dx 




dS 




'dt 



^dx J \9x 
-S{p++p-). 

We are interested in travelling wave solutions that develop from a pointwise 
inoculation of cells into a domain containing uniformly distributed nutrient 5*. 
In this scenario, p^ (defined as in each system) should form travelling pulses 
while S forms a travelling front and relevant boundary conditions will be 

± dp^ dS ^ 
Ox Ox 

S I as X ^ +00 , \ ■ ) 



S ^ S- as 



X — > — oo . 



Note that S- is currently unknown; we determine its value in the travelling 
wave analysis of Sectional Since p^ and S are physical quantities, we search 
for nonnegative travelling wave solutions, i.e. 

p^>0, S>0. 

It is clear that a travelling wave of this form cannot exist for S'c > 1 (extinction 
of bacteria) or for S'c < (infinite growth) and wc will therefore only consider 
systems that satisfy Sc € (0, 1). In the next section we analyse System (B) 
with respect to travelling wave solutions in order to obtain further insight. To 
do that, we use the rescaled system (|3.8p . 



4 Travelling wave analysis 

In this section we first apply the standard travelling wave ansatz to system 
(|3.8p and derive a necessary condition for the existence of non-negative trav- 
elling wave solutions. We then reduce the resulting ODE system to two com- 
ponents through a change of variables and utilizing an invariant manifold 
identified for the problem. Finally we use phase plane methods to analyse the 
existence and properties of travelling wave solutions. 



4.1 A necessary condition for the existence of travelling wave solutions 

Let us apply the travelling wave ansatz p^{x,t) = P^iO ~ p^{x — ct) and 
S{x,t) = S{^) = S{x — ct), where c is the unknown wave speed [29j . System 
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p.8p becomes 

(1 - c)(p+)' = - Ao (1 - X P+ + Ao (1 + X S') p- + (S- S,)p+ , 
-{l + c){p-y= Xoil-xS')p+ ~Xoil + xS')p- + {S-S,)p- , (4.1) 
-cS' ^-Sip+ +p-), 

where the primes denote derivatives with respect to the travelhng wave variable 
^. Note that any point on the S'-axis is a steady state of the system (|4.ip and 
that linear stability of such a steady state, {p'^,p~ ,S) = (0, 0, 5*,), is governed 
by the eigenvalues of the matrix A~^B, where 



'1-c 
A= I -1-c 1, B 
-cj 

The eigenvalues of A^^B are 



Ao + - 


- Sc Ao 





Ao 


— Ao + S", — 


ScO 










Ml = 0, ^i2,3 = — n — , (4.2) 

1 — 

where 

Ai{S,) = c'Xl + {S, -Sc~ 2Xo){S, - Sc). (4.3) 

Under the boundary conditions (|3.9p we look for nonnegativc solutions to (|4.1[) 
connecting steady states {p~^,p~,S) — (0, 0,5_) and {p~^,p~,S) = (0,0,1). To 
admit such a solution the latter must be a stable node, since a stable spiral 
would imply negative values for p^ . Hence, a necessary condition is Z\i (1) > 0, 
which is equivalent to 

c> c* ^ ^^{2Xo - I + Sc){l - Sc) . (4.4) 
Ao 

Given 2Ao > (1 — Sc) it is easy to show that c* <E [0, 1]. 

Theorem 1 A necessary condition for the existence of nonnegative travelling 
wave solutions of the system (j3.8p is 

2Xo > (1 - 5,). (4.5) 



The above condition is reasonable, as we expect the run duration to occur on 
a much faster time scale than proliferation processes. 
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4.2 Dimension reduction 

Let us now perform a change of variables by introducing the ccU density n — 
p+ + p~ and the cell flux j = ^ P - The travelling wave system (|4.ip can 
then be written as 

~ m' + f = {S - Sc)n , (4.6) 
-cj' + n' ^ 2AoX S'n +{S-Sc- 2Xo)j , (4.7) 
-cS" = -Sn , (4.8) 

where the boundary conditions for this system are 

. dn dj dS ^ ^ , 

",J,-^,Tp,-5 ^0 as ^->±oo, 

ox ox ox 

5 1 as ^ +00 , 

S ^ S- as ^ ^ -oo . 

From (|4.8p . we have Sn = cS' and, hence, n ~ c(ln S)'. Substituting into (|4.6p 
we obtain 

-cn' +j' = cS' -cScilnSy . 

Integrating and applying the boundary conditions at ^ — > +oo, an invariant 
manifold of the problem is given by 

—cn + j = c{S — 1) — cSc In S . 

With the definition f{S) = 5 — 1 — In S*, we obtain j = cn + cf{S), which 
can be used to eliminate j from the system (|4.6p - (|4.8p . For c 7^ 1 we can solve 
for n' and obtain the reduced system 

^^2^^ + 2 n{S -S,- Ao) + iS-Sc- 2Xo) f{S) 



(4.9) 



S'^-Sn. (4.10) 

c 



For c = 1 , we obtain 



Ao - + - v/(Ao -S + 5e)2 - 2Aox Sis - 5, - 2Ao)/(5) 



(4.11) 



2Aox^ 

S' = Sn, (4.12) 

where we chose the solution to the quadratic equation for n that satisfies the 
boundary conditions n — >■ as ^ — >■ ±00. 

It can be easily shown that f{S) = has two solutions in the region (0, 1] 
for all Sc e (0, 1) as follows. Since f'{S) = 1 — Sc/S, f{S) is monotonically 
decreasing for S G (0,5c) and monotonically increasing for S £ {Sc, 1]. With 
/(I) = 0, this implies f{Sc) < and, using f{S) — >■ 00 for 5 0, we obtain the 
existence and uniqueness of the second root of f{S) = 0: we call it € (0, Sc)- 
The existence of 5*1 and the negativity of f{S) for S G {Si, 1), together with the 
condition 2Ao > 1 — Sc, implies that n as given in (|4.1ip is positive everywhere, 
and that the given solution therefore satisfies the nonnegativity condition. 
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4.3 Steady states and their linear stability 

Using the two roots of f{S) — and under the condition (j4.5D . it is clear 
that there are two steady states of the system (|4.9|) - (j4.10p : {n, S) = (0, 1) and 
{n,S) = (OjS'i). Linearising the system (|4.9p - (|4.10p about its steady states 
generates a system of the form 

where, for the general steady state 5* € {^i, 1}, we have 



A = 



1 C L C O:^: 

^ 
\ c 

with 



2c 1 

traced = ^ (5** - S'c - Aq) , det^ = ^(5, - S'c - 2Ao)(S'* - S'c) 

1 — c^ 1 — c^ 



The eigenvalues of A are identical to ^2,3 as given in (|4.2p . The steady state 
(0, 1) is therefore a stable node for all c S (c*, 1) with c* as defined in (|4.4p . 
Similarly, it can be seen that the steady steady (0, Si) is a saddle point. The 
eigenvectors corresponding to the eigenvalues fi2.3 take the form 



S^ 



Vl,2 = [ ^2,3 ^ 

In the n — S plane, the slopes of the eigenvectors are given by 

fci,2(c) = ^. (4.13) 
For the steady state (n, 5*) = (0, 1) this slope ean be written in the form 

where we define A — c^Aq + (1 — Sc — 2Ao)(l — Sc) similarly to 



4.4 Case I: No chemotaxis (k = oo) 

We first consider the case where the chemotactic sensitivity x (given by (|3.6|) ) 
vanishes, i.e cells do not respond chemotactically to changes in S. Here, trav- 
elling waves are generated solely through proliferation of bacteria at the wave 
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Fig. 3 Illustration of the travelling wave solution calculated using the ODE system 
ll4:9ll- ll4J0ll for X = 0, Ao = 10, c = c* = 0.3122 and Sc = 0.5. 

(a) Trajectory of travelling wave solution. Solid line: trajectory, dashed line: n-nullcline, 
dotted line: circumference of invariant region Q introduced in the proof of Theorem [gl 

(b) Travelling wave solution in Solid line: n, dashed line: S. 



front. To understand the wave behaviour we perform a phase plane analysis 
for the ODE system (j4.9p - (|4.10p . Using k = oo (i.e. x — 0)i i* reduces to 



1 - c 

S' = —Sn . 
c 



2n{S -S,- Ao) + (5 - 5e - 2Ao)/(^) 



(4.15) 



Thus, the slope of a trajectory in the n — S plane can be written as 

dn _ 2n{S-Sc-Xo) + {S-Sc-2Xo)f{S) 

dS* ^ 1 - c2 Sn 

Additionally, an expression for the n— nuUcline Fn is given by 

_ S — Sc — 2Ao 
"""2(^-^,-Ao)^^'^^' 

and the S'-nullcline is simply 

n = 0, or 5 = 0. 

Let us now show that travelling waves exist for the reduced system (|4.15p . 

Theorem 2 For the case X = (which is equivalent to k = oo), a unique 
travelling wave solution for the system p.Sp exists for all c ^ (c*, 1). 

Proof For any c G (c*, 1) we can define a region fl (see Figure [3l[a)), enclosed 
by the line n ~ k2{S — 1) (with ^2 defined in (|4.14[) '). the S'-nullcline n — Q 
and the line S = Si. We will first show that f] is an invariant region of the 
system p.8p . Since S is non-decreasing everywhere in J7 and n' is non- negative 
for n = and S E [Si,l], we need only to show that the direction field on 
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the segment Fi = {{n,S) : n — k2{S — £ [5*1,1)} points from the top 
half of the plane above this segment towards the bottom. Since S is strictly 
increasing we require 

dn 



dS 



Indeed, 



< fc2 (< 0) , 



dn 
dS 



= 2- 



S-S, 



= 2 



S 

S — Sc — Xq 



\o^iS-Sc- 2Ao)/(5) 



< 2- 



S 

Sc 



Ao 



S 



S{S-l)k2 

{S-Sc- 2Ao)/(5) 
{S-Sc-2Xa){l- 



(1 



Sc^Xl 



Ao + VZ) ^ 



where we used (|4.14p in the first step and the relation f{S)/{S — 1) < 1 ^ •S'c 
for all S £ [5*1, 1]. Using the fact that k2 and {S ~ Sc~ 2Ao) are negative, we 
can use the definition of c* and the fact that 5 < 1 to obtain 



1 — dn 
c2 dS 



^ -2Ag + 3Ao(^ - Sc) -{S- Sc){l - Sc) 
5(2 Ao + - 1) 

2Ao ~l~ Sc — S 



A. 



< 



S{2Xo + Sc-l) 
-2Ag + 2Ao(5 - Sc) + Ao(l - Sc) - {S - 5c)(l - Sc) 
S{2Xo + Sc-l) 

2Ao ~\~ Sc — S 



2An + Sc-1' ' 

(2Ao + 5c-l)(Ao + 5c-5) _ 2Ao + Sc-S 
S{2Xo + Sc-l) 2Ao + Sc-l 



A, 



< -Ao + 1 - Sc~ VA = 



where we used S < 1 throughout the derivation. We can therefore conclude 
that i7 is an invariant region of the system p.Sp . Noting that at the steady state 
(n, 5*) = (0, 5i) the unstable manifold has a positive slope (fci.2 = M2.3C/S'*), 
i.e. it points into the region ]7, and using the fact that S is strictly increasing 
inside !7 for n > we can conclude that, for each c > c*, there is a heteroclinic 
orbit starting from (0, 5i) and finishing at (0, 1), corresponding to a travelling 
wave solution of the PDE system p.8|) . □ 



16 



Franz, Xuc, Painter, Erban 



4.5 Case II: Increasing chemotaxis (0 < k < oo) 

Decreasing n corresponds to an increase in the chemotactic sensitivity x ii^ 
the ODE system (j4.9[) - (j4.10[) and the slope of trajectories in the n — S plane 
is determined by 

dn c2 2n{S- Sc-\o) + {S- Sc -2\o)f{S) 2Aox 

zzz -4- 77 

dS* 1 - c2 Sn 1 - c2 ■ 

It is noted that the above slope is larger than that for the non-chemotaxis case 
within the region of interest n > 0. Due to this increase the region Q for the 
proof of Theorem 1 is no longer invariant for this system and a travelling wave 
solution to p.8p does not necessarily exist for all c g (c*, 1). The n-nuUcline 
for the full ODE system (j4.9[) - (|4.10p is given as the solution of the quadratic 
equation 

^^^n^ + 2c{S -Sc- Ao)?7. + c{S -Sc- 2\o) f{S) = 0. 
c 

For a given wave speed c, the rt-nuUcline can therefore be calculated as 



c(Ao + Sc-S)± ^/MS)\ , 
with 

A2{S) = c2(Ao + Sc- - 2Aox S{S - - 2\o) f{S) . 

We can see that A2{S) — > — oo as 5 oo due to its leading order term 
— 2Aox5''^. Therefore, as 5* becomes large, no n-nullcline exists and n' is posi- 
tive everywhere. Additionally, ^2(8) might have further roots and, in partic- 
ular, A2{S) might be negative in parts (or the whole) of region S € [5*1, 1]. 
This again means that n is strictly growing in these parts of the domain. 

We detect three different types of behaviours of trajectories starting close 
to {n, S) = (0, 5*1), plotted in Figure H) In particular, we can see each of these 
behavioural types for different values of x and despite different configurations 
of the nuUclines. In the top two plots of Figure H] we present the case of a 
diverging solution. Examining ODE (|4.9p . we observe that for large n, n grows 
quicker than 0{n^) and the divergence can be identified as a finite-time blow- 
up. In the second case, depicted in the two plots in the middle of Figure SI 
the trajectory converges to the steady state (0, 1), but does so after entering 
the region 5 > 1 and thereafter the region n < 0. Note that the steady state 
(0, 1) is still a stable node in this case and that this overshoot is therefore not 
a spiralling effect. Since these trajectories do not correspond to a non-negative 
solution of the ODE system (|4.9p - (|4.10p . they do not represent travelling wave 
solutions to the original problem. The last case, presented in the plots on the 
bottom of Figure m corresponds to an acceptable solution and is characterised 
by the convergence to (0, 1) without crossing the line S = 1. 
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X = l,c = 0.5884 x = 0.3, c = 0.3 




s s 



Fig. 4 Trajectories of the ODE system 114. 9l l - 114.1011 that highlight the three different cases. 
Parameters in all plots are Aq = 10, Sc = 0.5. Solid line: trajectory, dashed line: n-nullcline, 
dotted lines: n = and S = I. 



4.6 Case III: Infinite clieniotactic sensitivity (k ~ 0) 

As K decreases further we observe that the minimal wave speed necessary to 
allow a non-negative travelling wave solution of (|3.8p increases. In the limit 
K — > 0, the ODE system (j4.9p - (j4.10p no longer has convergent solutions. How- 
ever, in this limit the linearisation assumption leading to these ODEs and the 
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system (j3.8p is no longer valid and we must consider the original turning kernel 
as defined in p. 21) . In the limit k — > the turning rate in the hybrid model 
therefore becomes 

^-\2Ao, foryW<0. ^^-^^^ 

Hence, bacteria moving in a favourable direction do not turn, indicating that 
the wave speed achieved in this limit should evolve to c = s = 1. In [40] it was 
shown, for a slightly different turning kernel, that travelling waves can exist 
even without growth terms and that their wave speed satisfies c = s. 



5 Computational analysis of the wave speed 

In this section wc computationally compare wave speeds from the hybrid 
model with those of the fully continuous models. Specifically, wc investigate 
the regimes in which the latter provide an acceptable insight into the travel- 
ling wave behaviour of the hybrid model, and where they differ. We begin by 
investigating the non-chemotaxis case, where the minimum wave speed c* for 
the continuum systems was determined in (|4.4p . In Section 15.21 we show how 
the wave speed depends on the value of k, and correspondingly the chemotac- 
tic sensitivity x hi the macroscopic model. A comparison with hybrid models 
without cell proliferation is given in Section [^31 We conclude this section with 
a discussion into the effect and origin of oscillations observed under increasing 
the adaptation time ta- 



5.1 Case I: No chemotaxis (k oo) 

In Section ini wc analysed the macroscopic PDEs in the absence of chemotaxis. 
Travelling wave solutions were shown to exist for all wave speeds c S (c*, 1), 
with c* determined by (|4.4p . In Figure [5I^a), variation of (|4.4p as a function of 
Ao is illustrated; we note that wave speeds determined through simulation of 
the PDE systems correspond exactly (to accuracy of the numerical approxi- 
mation) with the analytical wave speeds. We now numerically investigate the 
wave speed for the case x = in the hybrid model. 

For our simulations we consider the same parameters and methods as de- 
scribed in Section 12.31 specifically, we set the system parameters Sc ~ 0.5, 
s = 1 and Ds = 0. For the computations we consider a time step At = 10~^, 
a spatial resolution of Ax = 0.25 on a domain with length L = 100, and sim- 
ulate the system until the value of S" at x = 60 falls below 0.5. The profiles at 
this time, together with the time when 5 at a; = 20 falls below 0.5, are used 
to estimate the wave speed. 

The measured wave speed for varying Aq is illustrated in Figure [Slja), along 
with c* as predicted from the travelling wave analysis. While the relationship 
is similar in shape, wc note that at all values of Aq tested the measured wave 
speed lies below the analytical value c* . In the literature it has been observed 
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Fig. 5 Measured wave speed in the hybrid model. Crosses: individual simulations, 
ensemble averages. Parameters are as described in the text. 

(a) Wave speed in dependence of Xq for Nq = 10,000. Dashed line: c* given by 1 14.411 . 

(b) Wave speed as a function of No with Ao = 10. Dashed line: c* computed by 114.411 . 



dots: 



that inaccuracies in numerical schemes can lead to an increase in wave speeds 
[33| , therefore rendering the lower wave speed seen in Figure [SJa) as counter 
intuitive. 

Nevertheless, we can provide the following explanation for the distinct val- 
ues in the continuum and hybrid models. For the zero-chemotaxis case, wave 
generation and movement is solely determined by growth ahead and death 
behind the wave. In the continuum model an outermost "fractional bacteria 
population" can extend significantly beyond the wave front, since some pro- 
portion of the initial population never turns left, and hence far into the region 
where S is very close to its initial value of 1. Yet this fractional population still 
grows exponentially {dp^/dt « {1 — Sc)p^), seeding the growth and expansion 
of the population. The finite/discrete nature of the hybrid model precludes any 
fractional bacterium: the forward "tail" is necessarily finite and growth will 
not occur beyond the outermost individual. 

For the above explanation to hold we would expect a dependence of the 
measured wave speed on the initial number of bacteria A^o^ continuous densi- 
ties provide a closer approximation under larger numbers of bacteria and we 
would expect convergence in the wave speed to c*. Simulations in Figure [SJb) 
demonstrate this property, corroborating our interpretation. 



5.2 Case II: Increasing chemotaxis (0 < k < oo) 

In the second set of numerical experiments we measure the dependency of 
the wave speed on the critical parameter k, i.e. we determine the effect of 
increasing chemotaxis as k decreases. We compare the results measured for 
the hybrid system with the continuous Systems (A) and (B). 

We use the same parameters as in Section 15.11 and results are shown in 
Figure [51 The results demonstrate the regimes where correspondence across 
the varying modelling levels occurs: while the hybrid model (dotted line) cor- 
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Fig. 6 Comparison between wave speeds of the various models in dependence of k. Dot- 
ted line: hybrid model, red solid line: mesoscopic System (A), dashed line: linearised Sys- 
tem (B). Parameters are as described in the text. 

responds well with its closest continuous version (mesoscopic System (A), red 
solid line) over a wide range of k, it only corresponds with System (B) (black 
dashed line) for larger k, diverging as k decreases. Note that the turning rate 
(|3.6p used for System (B) becomes negative at small values of k and we limit 
the range of k studied accordingly. 

At larger k all three models converge to a value close to c* as k grows: in 
this regime the main assumption proposed for the linearisation (|S'(a;) — y^^^ | <C 
n) holds and we obtain good quantitative agreement. While this assumption 
becomes less acceptable as we decrease k, leading to the divergent behaviour 
described above, we note that all models show the same qualitative agreement: 
increasing chemotactic responses leads to an increase in the wave speed. Note 
that the results for System (B) can be identically replicated using the ODE 
system (|4.9p - (|4.10[) and a search algorithm for the smallest value of c that 
admits a nonnegative solution to the system. 

These numerical experiments demonstrate that chemotaxis has a signifi- 
cant effect on the speed of movement and that the waves cannot solely be 
explained by growth and death terms. Rather, we interpret birth and death 
processes as stabilisers to what would otherwise be transient waves [T5ll40j . 
This interpretation is in agreement with the results presented in Figure [H as 
the initial wave speed for the system without growth seems to be similar to 
the wave speed of the system including growth and death terms. 



5.3 Oscillations in the wave speed 



An additional observation wc made during the numerical experiments of the 
hybrid model is that for increasing values of the adaptation time ia, the wave 



Travelling waves in hybrid chemotaxis models 



21 



speed starts to differ strongly from the mesoscopic System (A) , an effect that 
we identified to be due to oscillations in the behaviour of the wave. In Fig- 
ure El^a) we present an example of strongly oscillating wave speeds (where the 
wave speed is measured as rate of change of the average position of bacteria). 
This example occurred for the parameters Sc = 0.5, Aq = 10, k = 0.001 and 
ta ~ 4:. We can also clearly see that the wave speed is correlated to the current 
number of agents in the system. In the literature similar effects of oscillating 
waves in stochastic models have been observed [28ll32] . 

In Figure [TJb), we present the form of the wave at different times through- 
out the simulation. It is clearly visible that the shape differs significantly at 
different times. One reason these oscillations occur when ta is very high is 
that a bacterium that happens to be in front of the wave experiences a very 
high value of S, whilst its internal dynamics only adapt very slowly. This, 
in combination with the low value of k, leads to a bacterium that does not 
switch direction for a long time and will proliferate at a high rate. This implies 
that a spike of bacteria forms in front of the wave that moves faster than the 
rest of the wave. We can clearly see such a spike in the left-most waveform 
in Figure [71[b) . Once the frontrunning bacterium and its copies have turned, 
the wave goes into a reordering phase (second and third waveform), until, 
eventually, a new spike emerges (4th waveform). 

In Figure EJ^c) we plot the wave speed over time for a smaller value of 
ta- We can see that the oscillations are less severe and more frequent than in 
Figure [71[a), which is in agreement with the explanation above. As we decrease 
ta the frontrunning bacteria will adapt quicker to their surroundings and are 
thereby more likely to turn. We show the influence of changing Nq on the 
oscillating behaviour in Figure El^d) . The oscillations seem to occur with a 
similar frequency but more regular to those before, which can be explained 
by the increased likelihood of frontrunning bacteria with a higher number of 
agents and reduced noise in the system. 



6 Discussion 

In this paper we presented a hybrid model of chemotaxis, incorporating a bio- 
logically realistic turning kernel introduced in |3D] . We analysed the travelling 
wave behaviour of this hybrid system using mesoscopic and macroscopic equa- 
tions, deriving an analytical value for the expected wave speed in the case 
of no chemotaxis. As chemotaxis increases we demonstrated (analytically and 
numerically) that the expected wave speed increases, indicating that the wave 
that forms is not solely driven by growth and death processes. In contrast to 
the transient waves observed for the hybrid model in the absence of growth 
and death terms |15j . the (numerical) waves observed here in their presence 
are stable, indicating the stabilising effect of birth and death. The numerical 
analysis reveals that the macroscopic equations derived through linearisation 
of the turning kernel can qualitatively describe the change in wave speed as 
chemotaxis increases, but that there are significant quantitative differences 
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Fig. 7 Oscillations in the wave speed of the hybrid model | |2. 21 1 - 112.61 1 and 112.1011 . 

(a) Wave speed in comparison to current number of particles for No = 10, 000, ta = 4. Solid 
line: wave speed, dashed line: number of particles, dotted lines: times of wave forms shown 
in panel (b). 

(b) Waveform at 4 distinct times marked in panel (a) from left to right. 

(c) As in (a) with Nq = 10, 000, ta = 2. 

(d) As in (a) with No = 50, 000, ta = 4. 
Other parameters are given in Section 15.31 



between the two systems. Additionally, we observed oscillations in the wave 
movement, an effect that had been seen in similar systems in the literature 
[52] and that cannot be explained using mean-field approximations. 

To date, travelling waves in chemotaxis models have mainly been analysed 
from the perspective of macroscopic PDE models of chemotaxis [TMT5] . The 
existence of travelling waves for continuum models with growth terms is well 
established [36. 30.1 22] . While hybrid models have been used to study pattern 
formation in bacterial chemotaxis [T71I55] , these studies have not analysed the 
travelling wave patterns observed in bacterial cell populations. 

Recently, experimental studies using microfluidic techniques tracked cell 
trajectories within a traveling pulse, and revealed that persistence of direction 
in cell movement accounts for 30% of the macroscopic speed of the travel- 
ing pulse |35| . The hybrid model framework studied here provides a natural 
method for direct comparison of model predictions with experimental mea- 
surements of cell trajectory, and this is left as future work. 
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